23. Ordinal Predicted Variable

23.1 Modelling Ordinal Data with an Underlying Metric Variable

• Ordinal data is often a bracketing/bucketing of an underlying metric value (which may be normally distributed)
• Even if the underlying data is normally distributed the ordinal values may not resemble a normal distribution (Fig 23.1)
• The probability of a particular ordinal outcome is the area under the normal curve between the thresholds of that outcome.
  • Can be calculated using differences in the cumulative normal distribution: \(p(y=k|\mu,\sigma,\{ \theta_j \})=\Phi((\theta_k-\mu)/\sigma)-\Phi((\theta_{k-1}-\mu)/\sigma)\) where \(\theta_k\) are the thresholds of the outcomes.
• Takeaway: we might be able to model a distribution of ordinal values as consecutive intervals of a normal distribution on an underlying metric scale with appropriately partitioned thresholds.
  • BB: this doesn’t really help you in modelling if the thresholds can be chosen freely? If they cannot be chosen freely then we are applying some strange constraints and we usually would not really have much information on how to apply these or why they would be a good model of the data?

23.2 The Case of a Single Group

• If there are K ordinal values, the model has K+1 parameters: \(\theta_1,...,\theta_{K-1},\mu,\sigma\).
  • If we increase all the thresholds by \(a\) that gives the same result as decreasing \(\mu\) by \(a\)
  • Similar with expanding or contracting the thresholds and varying \(\sigma\)
  • No uniquely correct choice of how to “pin down” the axis by setting two of the parameters, we choose:
    • \(\theta_1 = 1 + 0.5\)
    • \(\theta_{K-1} = K - 0.5\)
    • BB: end up with mean and standard deviation effectively being controlled by amount of data in the edge groups, then the middle groups fit themselves around what’s left..?

"
model {
  for (i in 1:Ntotal) {
    y[i] ~ dcat(pr[i,1:nYlevels])
    pr[i,1] <- pnorm(thresh[1], mu, 1/sigma^2)
    for (k in 2:(nYlevels-1)) {
      pr[i,k] <- max(0, pnorm(thresh[k], mu, 1/sigma^2) - pnorm(thresh[k-1], mu, 1/sigma^2))
    }
    pr[i,nYlevels] <- 1 - pnorm(thresh[nYlevels-1], mu, 1/sigma^2)
  }
}
"

"
mu ~ dnorm((1_nYlevels)/2, 1/(nYlevels)^2)
sigma ~ dunif(nYlevels/1000, nYlevels*10)
for (k in 2:(nYlevels-2)) {
  thresh[k] ~ dnorm(k+0.5, 1/2^2)
}
"

23.3 The Case of Two Groups

Example: ask about Happiness for two groups - one after watching sports car ads, one after watching charity ads. Assume that both groups access the same underlying metric scale and thresholds, what differs is the mean and variance of the feeling.

23.4 The Case of Metric Predictors

Figure 23.6 for important explanation of how we do linear regression to predict the underlying metric variable, then we apply our same model to convert underlying metric variable to ordinal predicted variable. Ordinal probit regression.

Figures 23.7 to 23.10 interesting examples.

23.5 Posterior Prediction

How are posterior predicted probabilities with 95% HDI computed? At each step in MCMC chain compute predicted outcome probabilities using equations given and explained. From the full chain of credible posterior predicted outcome probabilities the central tendency and 95% HDI are then computed. Brief example how to do this below: use the MCMC chain coda object returned from JAGS (codaSample) then do the processing of it in R.

mcmcMat = as.matrix(codaSamples)
chainLength = nrow(mcmcMat)
beta0 = mcmcMat[,"beta0"]
beta1 = mcmcMat[,"beta"]
thresh = mcmcMat[,grep("^thresh",colnames(mcmcMat))]
sigma = mcmcMat[,"sigma"]

outProb = matrix(NA,nrow=chainLength,ncol=max)
for (stepIdx in 1:chainLength) {
  mu = beta0[stepIdx] + beta1[stepIdx] * xProbe
  threshCumProb = pnorm((thres[stepIdx, ] - mu)/sigma[stepIdx])
  outProb[stepIdx, ] = c(threshCumProb,1) - c(0,threshCumProb)
}
outMed = apply(outProb, 2, median, na.rm=True)
outHdi = apply(outProb, 2, HDIofMCMC)

23.6 Generalizations and Extensions

Approaches to outliers:

• Use heavy-tailed distribution instead of normal distribution to describe noise. Can use cumulative t function which is build into JAGS as pt (analogous to pnorm).
• Model the data as being from two different sources as we did in Robust Logistic Regression; one source is our cumulative normal model and the other is “guessing”; so we mix our models predicted probabilities with some random guessing probability.

Variable selection: just as predictors in linear regression or logistic regression can be given inclusions parameters, so can our predictors in thresholded cumulative-normal regression.

Nominal predictors: works just fine, look back to Figure 21.12 to see how this might work. The only change is putting thresholds on the normal noise distribution to create probabilities for the ordinal outcomes.

23.7 Exercises

setwd("./DBDA2Eprograms")
source("DBDA2E-utilities.R") # Load definitions of graphics functions etc.
source("Jags-Yord-XmetMulti-Mnormal-Guessing-Example.R")

setwd("./DBDA2Eprograms")
source("DBDA2E-utilities.R") # Load definitions of graphics functions etc.
source("Jags-Yord-XmetMulti-Mnormal-tDist-Example.R")